***************************************************************
* PROGRAM TO RUN AN INDIVIDUAL SET OF 2SLS REGRESSION RESULTS *
***************************************************************

capture program drop directctrl2sls
program define directctrl2sls

	set more off

	syntax [if/], outcome(varname) [ctrlvars(namelist)] trimlevel(numlist max=1 >=0 <100) [header fh(name) printtoscreen]
	
	* Create list of control variables	
	local ctrlfvarlist ""
	foreach var of local ctrlvars {
		forval lag = 1/`lags' {
			if "`var'" == "survives" {
				local ctrlfvarlist "`ctrlfvarlist' i.`var'_l`lag'"
			}
			else local ctrlfvarlist "`ctrlfvarlist' c.`var'_l`lag'"
		}
	}

	if (("`ctrlvars'" == "all") | (("`ctrlvars'" == "survives") & ("`outcome'" == "survives"))) local lastpre = 20 - `lags' - 1
	else local lastpre = 20 - 1

	* Work out if condition
	if "`if'" == "" local if "if 1"
	else local if "if (`if')"
	local if "`if' & (audityear >= 1999+4) & (yrsfiling >= 2)"

	* Work out trimming
	if (`trimlevel' == 0 ) {
		gen byte trimmed = 0
	}
	else {	//if trimlevel>0, construct variable noting which obs are trimmed
		* Cut (exclude) above a given percentile
		* Calculate percentile to trim at
		qui gquantiles `outcome' if (`outcome' > 0 & `outcome' < .), _pctile percentiles(`=100 - `trimlevel'')
		
		* Create trimmed variable
		gegen byte trimmed = max(`outcome' > `r(r1)'), by(utr_no audityear)
		
	}
	
	* Modify if condition to avoid having the same variable on both sides of the equals sign
	if ((inlist("`ctrlvars'","all","own")) | (("`ctrlvars'" == "survives") & ("`outcome'" == "survives"))) local if "`if' & (!inrange(yrssince,-4,-1))"
	* Modify if condition for trimming
	local if "`if' & (trimmed == 0)"

	* Create differenced outcome variable
	gegen double `outcome'_avg4 = total(`outcome'*inrange(yrssince,-4,-1)/4), by(utr_no audityear)
	qui replace `outcome'_avg4 = . if minyrssince > -4
	qui gen double `outcome'_diff4 = `outcome' - `outcome'_avg4
	drop `outcome'_avg4
	local outcomename "`outcome'_diff4"
	
	* 1st stage regression
	**********************
	
	* Set up the matrices
	matrix stage1b = J(1,21,.)
	matrix stage1V = J(1,21,.)
	* Label the column names
	local colnames ""
	forvalues yr20 = 10/30 {
		local colnames "`colnames' `yr20'.yrssince20#1.treatment"
	}
	matrix colnames stage1b = `colnames'
	matrix colnames stage1V = `colnames'
	
	* Run the 1st stage regression
	foreach yr20 of numlist 10/`lastpre' 20/30 {
		gen byte endogvar = i`yr20'.yrssince20*actuallyaudited
		qui reg endogvar i(10/`lastpre' 20/30).yrssince20 i(1999/2012).tax_year `ctrlfvarlist' i(10/`lastpre' 20/30).yrssince20#i.treatment `if', vce(cluster utr_no)
		matrix eb = e(b)
		matrix eV = e(V)
		local ebcolnumb = colnumb(eb,"`yr20'.yrssince20#1.treatment")
		matrix stage1b[1,`yr20'-10+1] = eb[1, `ebcolnumb']
		matrix stage1V[1,`yr20'-10+1] = eV[`ebcolnumb', `ebcolnumb']
		drop endogvar
		matrix drop eb eV
	}
	
	* 2SLS regression
	di `"qui ivregress 2sls `outcomename' i(10/`lastpre' 20/30).yrssince20 i(1999/2012).tax_year `ctrlfvarlist' (i(10/`lastpre' 20/30).yrssince20#i.actuallyaudited=i(10/`lastpre' 20/30).yrssince20#i.treatment) `if', vce(cluster utr_no)"'
	qui ivregress 2sls `outcomename' i(10/`lastpre' 20/30).yrssince20 i(1999/2012).tax_year `ctrlfvarlist'  (i(10/`lastpre' 20/30).yrssince20#i.actuallyaudited=i(10/`lastpre' 20/30).yrssince20#i.treatment) `if', vce(cluster utr_no)
	
	* Store results
	matrix eb = e(b)
	matrix eV = e(V)

	*make it print nicely
	if (("`ctrlvars'" == "") local ctrlvars "NA"
	
	* Print header if requested
	if ("`header'" == "header") {
		if ("`printtoscreen'" == "printtoscreen") di as text %12s "method" %24s "outcome" %12s "ctrlvars" %12s "trimlevel" %10s "yrssince"%14s "1ststage" %14s "1ststagese" %14s "1ststagepval" %14s "treateff" %14s "treateffse" %14s "pvalue" %14s "N"
		if ("`fh'" != "") file write `fh' "method,outcome,ctrlvars,trimlevel,yrssince,1ststage,1ststagese,1ststagepval,treateff,treateffse,pvalue,N" _n
	}
	

	* Print results
	forvalues yr20 = 15/28 {
		local ebcolnumb = colnumb(eb,"`yr20'.yrssince20#1.actuallyaudited")
		if ("`printtoscreen'" == "printtoscreen") {
			# delimit ;
			di as text %12s "late" %24s "`outcome'" %12s "`ctrlvars'" %12.0g `trimlevel' %10.0g `yr20'-20
				%14.3f stage1b[1, `yr20'-10+1] %14.3f sqrt(stage1V[1,`yr20'-10+1])
				%14.3f 2*ttail(e(N_clust), abs(stage1b[1, `yr20'-10+1]/sqrt(stage1V[1,`yr20'-10+1])))
				%14.3f eb[1, `ebcolnumb'] %14.3f sqrt(eV[`ebcolnumb',`ebcolnumb'])
				%14.3f 2*ttail(e(N_clust), abs(eb[1, `ebcolnumb']/sqrt(eV[`ebcolnumb',`ebcolnumb'])))
				%14.0g e(N)
				;
			# delimit cr
		}
		# delimit ;
		if ("`fh'" != "") file write `fh' "late,`outcome',`ctrlvars'," (`trimlevel') "," (`yr20'-20) ","
			(stage1b[1, `yr20'-10+1]) "," (sqrt(stage1V[1,`yr20'-10+1])) ","
			(2*ttail(e(N_clust), abs(stage1b[1, `yr20'-10+1]/sqrt(stage1V[1,`yr20'-10+1])))) ","
			(eb[1, `ebcolnumb']) "," (sqrt(eV[`ebcolnumb',`ebcolnumb'])) ","
			(2*ttail(e(N_clust), abs(eb[1, `ebcolnumb']/sqrt(eV[`ebcolnumb',`ebcolnumb'])))) ","
			(e(N)) _n
			;
		# delimit cr
	}
	matrix drop eb eV
	
	drop trimmed

	drop `outcomename'

end
